#include "precompiled.h"
#pragma hdrstop

#include <cstdio>
#include <climits>

#if (INF_COMPILER==COMPILER_MSVC) && (INF_COMP_VER>=1300)
#define sprintf sprintf_s
#endif

namespace inf {
    namespace math {
        const float LIMITS_FLOAT_MAX        = 3.402823466e+38F;
        const float LIMITS_FLOAT_MIN        = 1.175494351e-38F;
        const double LIMITS_DOUBLE_MAX      = 1.7976931348623158e+308;
        const double LIMITS_DOUBLE_MIN      = 2.2250738585072014e-308;

        const int IEEE_FLT_MANTISSA_BITS	= 23;
        const int IEEE_FLT_EXPONENT_BITS	= 8;
        const int IEEE_FLT_EXPONENT_BIAS	= 127;
        const int IEEE_FLT_SIGN_BIT			= 31;

        const int IEEE_DBL_MANTISSA_BITS	= 52;
        const int IEEE_DBL_EXPONENT_BITS	= 11;
        const int IEEE_DBL_EXPONENT_BIAS	= 1023;
        const int IEEE_DBL_SIGN_BIT			= 63;

        const int IEEE_DBLE_MANTISSA_BITS	= 63;
        const int IEEE_DBLE_EXPONENT_BITS	= 15;
        const int IEEE_DBLE_EXPONENT_BIAS	= 0;
        const int IEEE_DBLE_SIGN_BIT		= 79;

        const float	math::PI				= 3.14159265358979323846f;
        const float	math::TWO_PI			= 2.0f * PI;
        const float	math::HALF_PI			= 0.5f * PI;
        const float	math::ONEFOURTH_PI		= 0.25f * PI;
        const float	math::THREEFOURTHS_PI	= 0.75f * PI;
        const float math::E					= 2.71828182845904523536f;
        const float math::SQRT_TWO			= 1.41421356237309504880f;
        const float math::SQRT_THREE		= 1.73205080756887729352f;
        const float	math::SQRT_1OVER2		= 0.70710678118654752440f;
        const float	math::SQRT_1OVER3		= 0.57735026918962576450f;
        const float	math::M_DEG2RAD			= PI / 180.0f;
        const float	math::M_RAD2DEG			= 180.0f / PI;
        const float	math::M_SEC2MS			= 1000.0f;
        const float	math::M_MS2SEC			= 0.001f;
        const float	math::INFINITY			= 1e30f;
        const float math::FLOAT_EPSILON		= 1.192092896e-07f;
#ifdef INF_WIN_X86_SSE
        const float math::SSE_FLOAT_ZERO	= 0.0f;
        const float math::SSE_FLOAT_255		= 255.0f;
#endif

        namespace Private {
            bool		initialized		= false;
            dword		iSqrt[SQRT_TABLE_SIZE];		// inverse square root lookup table
        }

        /*
        ===============
        math::init
        ===============
        */
        void math::init( void ) {
            union Private::_flint fi, fo;

            for ( int i = 0; i < Private::SQRT_TABLE_SIZE; i++ ) {
                fi.i	 = ((Private::EXP_BIAS-1) << Private::EXP_POS) | (i << Private::LOOKUP_POS);
                fo.f	 = (float)( 1.0 / sqrt( fi.f ) );
                Private::iSqrt[i] = ((dword)(((fo.i + (1<<(Private::SEED_POS-2))) >> Private::SEED_POS) & 0xFF))<<Private::SEED_POS;
            }

            Private::iSqrt[Private::SQRT_TABLE_SIZE / 2] = ((dword)(0xFF))<<(Private::SEED_POS); 

            Private::initialized = true;
        }

        /*
        ================
        math::float_to_bits
        ================
        */
        int math::float_to_bits( float f, int exponentBits, int mantissaBits ) {
            int i, sign, exponent, mantissa, value;

            assert( exponentBits >= 2 && exponentBits <= 8 );
            assert( mantissaBits >= 2 && mantissaBits <= 23 );

            int maxBits = ( 1 << ( exponentBits + mantissaBits ) ) - 1;
            int minBits = 1;

            float max = bits_to_float( maxBits, exponentBits, mantissaBits );
            float min = bits_to_float( minBits, exponentBits, mantissaBits );

            if ( f >= 0.0f ) {
                if ( f >= max ) {
                    return maxBits;
                } else if ( f == min ) {
                    return minBits;
                } else if ( f < min ) {
                    return 0;
                }
            } else {
                if ( f <= -max ) {
                    return ( maxBits | ( 1 << ( exponentBits + mantissaBits ) ) );
                } else if ( f == -min ) {
                    return ( minBits | ( 1 << ( exponentBits + mantissaBits ) ) );
                } else if ( f > -min ) {
                    return 0;
                }
            }

            i = *reinterpret_cast<int *>(&f);
            sign = ( i >> IEEE_FLT_SIGN_BIT ) & 1;
            exponent = ( ( i >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
            mantissa = i & ( ( 1 << IEEE_FLT_MANTISSA_BITS ) - 1 );
            value = sign << ( exponentBits + mantissaBits );
            value |= ( exponent + ( ( 1 << ( exponentBits - 1 ) ) - 1 ) ) << mantissaBits;
            value |= mantissa >> ( IEEE_FLT_MANTISSA_BITS - mantissaBits );
            return value;
        }

        /*
        ================
        math::bits_to_float
        ================
        */
        float math::bits_to_float( int i, int exponentBits, int mantissaBits ) {
            int sign, exponent, mantissa, value;

            assert( exponentBits >= 2 && exponentBits <= 8 );
            assert( mantissaBits >= 2 && mantissaBits <= 23 );

            if ( ( i & ( ( 1 << ( exponentBits + mantissaBits ) ) - 1 ) ) == 0 ) {
                // if exponent & mantissa are zero, value is plus or minus zero
                return 0.0f;
            }

            sign = i >> ( exponentBits + mantissaBits );
            exponent = ( ( i >> mantissaBits ) & ( ( 1 << exponentBits ) - 1 ) ) - ( ( 1 << ( exponentBits - 1 ) ) - 1 );
            mantissa = ( i & ( ( 1 << mantissaBits ) - 1 ) ) << ( IEEE_FLT_MANTISSA_BITS - mantissaBits );
            value = sign << IEEE_FLT_SIGN_BIT | ( exponent + IEEE_FLT_EXPONENT_BIAS ) << IEEE_FLT_MANTISSA_BITS | mantissa;
            return *reinterpret_cast<float *>(&value);
        }

        /*
        ================
        math::test_float_bit_conversions
        ================
        */
        void math::test_float_bit_conversions( void ) {
            // 
            // test every value that can be represented within the valid ranges of bits
            char testBuf[ 512 ];
            int numFails = 0;
            for ( int exponentBits = 2; exponentBits <= 8; exponentBits++ ) {
                for ( int mantissaBits = 2; mantissaBits <= 23; mantissaBits++ ) {

                    sprintf( testBuf, "Float Testing: exp %i mant %i ...... ", exponentBits, mantissaBits );
#if (CORE_COMPILER==COMPILER_MSVC)
                    OutputDebugString( testBuf );
#else
                    printf( "%s", testBuf );
#endif

                    int maxBits = ( 1 << ( exponentBits + mantissaBits ) ) - 1;
                    int minBits = 1;
                    float max = math::bits_to_float( maxBits, exponentBits, mantissaBits );
                    float min = math::bits_to_float( minBits, exponentBits, mantissaBits );

                    bool failed = false;

                    // test that all representable values convert to & from a float & maintain value
                    unsigned int totalNumbersRepresentable = 1 << ( exponentBits + mantissaBits + 1 );
                    for ( unsigned int i = 0; i < totalNumbersRepresentable; i++ ) {
                        // skip anything with zero mantissa & exponent
                        float testFloat = math::bits_to_float( *( ( int* )&i ), exponentBits, mantissaBits );
                        unsigned int testBits = math::float_to_bits( testFloat, exponentBits, mantissaBits );

#pragma warning( push )
#pragma warning( disable: 4389 )
                        if ( i == 1 << ( exponentBits + mantissaBits ) ) {
#pragma warning( pop )
                            // negative zero should be the same as positive zero (not IEEE-754, but simple)
                            if ( testBits != 0 ) {
                                failed = true;
                                break;
                            }
                        } else if ( testBits != i ) {
                            failed = true;
                            break;
                        }
                    }

                    for ( float i = -max * 2.0f; i < max * 2.0f; ) {
                        unsigned int testBits = math::float_to_bits( i, exponentBits, mantissaBits );
                        float testFloat = math::bits_to_float( testBits, exponentBits, mantissaBits );

                        // it will have already passed the test of bit->float for this value, if
                        // float_to_bits has returned anything in the valid range
                        unsigned int minForThisexp = testBits - ( testBits & mantissaBits ) + 1;
                        float minError = math::fabs( math::bits_to_float( minForThisexp, exponentBits, mantissaBits ) );

                        if ( i < -max ) {
                            if ( testFloat != -max ) {
                                failed = true;
                                break;
                            }
                        } else if ( i > max ) {
                            if ( testFloat != max ) {
                                failed = true;
                                break;
                            }
                        } else if ( i > -min && i < min ) {
                            if ( testFloat != 0.0f ) {
                                failed = true;
                                break;
                            }
                        } else {
                            float error = math::fabs( i - testFloat );
                            if ( error > minError ) {
                                failed = true;
                                break;
                            }
                        }

                        i += minError;
                    }

                    if ( failed ) {
                        numFails++;
                        sprintf( testBuf, "FAILED\n" );
                    } else {
                        sprintf( testBuf, "PASSED\n" );
                    }
#if (CORE_COMPILER==COMPILER_MSVC)
                    OutputDebugString( testBuf );
#else
                    printf( "%s", testBuf );
#endif
                }
            }
            sprintf( testBuf, "%i tests failed\n", numFails );

#if (CORE_COMPILER==COMPILER_MSVC)
            OutputDebugString( testBuf );
#else
            printf( "%s", testBuf );
#endif
        }

        float math::inv_sqrt( float x )
        {
            assert( x != 0.f );
#if 1
            return 1 / sqrtf( x );
#else
            dword a = ((union _flint*)(&x))->i;
            union _flint seed;

            assert( initialized );

            double y = x * 0.5f;
            seed.i = (( ( (3*EXP_BIAS-1) - ( (a >> EXP_POS) & 0xFF) ) >> 1)<<EXP_POS) | iSqrt[(a >> (EXP_POS-LOOKUP_BITS)) & LOOKUP_MASK];
            double r = seed.f;
            r = r * ( 1.5f - r * r * y );
            r = r * ( 1.5f - r * r * y );
            return (float) r;
#endif
        }

        double math::int_sqrt64( float x )
        {
            assert( x != 0.f );
#if 1
            return 1 / sqrtf( x );
#else
            dword a = ((union _flint*)(&x))->i;
            union _flint seed;

            assert( initialized );

            double y = x * 0.5f;
            seed.i = (( ( (3*EXP_BIAS-1) - ( (a >> EXP_POS) & 0xFF) ) >> 1)<<EXP_POS) | iSqrt[(a >> (EXP_POS-LOOKUP_BITS)) & LOOKUP_MASK];
            double r = seed.f;
            r = r * ( 1.5f - r * r * y );
            r = r * ( 1.5f - r * r * y );
            r = r * ( 1.5f - r * r * y );
            return r;
#endif
        }

        double math::sqrt64( float x )
        {
            return sqrt( x );
        }

        float math::sin( float a )
        {
            return sinf( a );
        }

        float math::sqrt( float x )
        {
            return sqrtf( x );
        }

        float math::sin16( float a )
        {
            float s;

            if ( ( a < 0.0f ) || ( a >= TWO_PI ) ) {
                a -= floorf( a / TWO_PI ) * TWO_PI;
            }
#if 1
            if ( a < PI ) {
                if ( a > HALF_PI ) {
                    a = PI - a;
                }
            } else {
                if ( a > PI + HALF_PI ) {
                    a = a - TWO_PI;
                } else {
                    a = PI - a;
                }
            }
#else
            a = PI - a;
            if ( fabs( a ) >= HALF_PI ) {
                a = ( ( a < 0.0f ) ? -PI : PI ) - a;
            }
#endif
            s = a * a;
            return a * ( ( ( ( ( -2.39e-08f * s + 2.7526e-06f ) * s - 1.98409e-04f ) * s + 8.3333315e-03f ) * s - 1.666666664e-01f ) * s + 1.0f );
        }

        double math::sin64( float a )
        {
            return sin( a );
        }

        float math::cos( float a )
        {
            return cosf( a );
        }

        float math::cos16( float a )
        {
            float s, d;

            if ( ( a < 0.0f ) || ( a >= TWO_PI ) ) {
                a -= floorf( a / TWO_PI ) * TWO_PI;
            }
#if 1
            if ( a < PI ) {
                if ( a > HALF_PI ) {
                    a = PI - a;
                    d = -1.0f;
                } else {
                    d = 1.0f;
                }
            } else {
                if ( a > PI + HALF_PI ) {
                    a = a - TWO_PI;
                    d = 1.0f;
                } else {
                    a = PI - a;
                    d = -1.0f;
                }
            }
#else
            a = PI - a;
            if ( fabs( a ) >= HALF_PI ) {
                a = ( ( a < 0.0f ) ? -PI : PI ) - a;
                d = 1.0f;
            } else {
                d = -1.0f;
            }
#endif
            s = a * a;
            return d * ( ( ( ( ( -2.605e-07f * s + 2.47609e-05f ) * s - 1.3888397e-03f ) * s + 4.16666418e-02f ) * s - 4.999999963e-01f ) * s + 1.0f );
        }

        double math::cos64( float a )
        {
            return cos( a );
        }

        void math::sincos16( float a, float &s, float &c )
        {
            float t, d;

            if ( ( a < 0.0f ) || ( a >= math::TWO_PI ) ) {
                a -= floorf( a / math::TWO_PI ) * math::TWO_PI;
            }
#if 1
            if ( a < PI ) {
                if ( a > HALF_PI ) {
                    a = PI - a;
                    d = -1.0f;
                } else {
                    d = 1.0f;
                }
            } else {
                if ( a > PI + HALF_PI ) {
                    a = a - TWO_PI;
                    d = 1.0f;
                } else {
                    a = PI - a;
                    d = -1.0f;
                }
            }
#else
            a = PI - a;
            if ( fabs( a ) >= HALF_PI ) {
                a = ( ( a < 0.0f ) ? -PI : PI ) - a;
                d = 1.0f;
            } else {
                d = -1.0f;
            }
#endif
            t = a * a;
            s = a * ( ( ( ( ( -2.39e-08f * t + 2.7526e-06f ) * t - 1.98409e-04f ) * t + 8.3333315e-03f ) * t - 1.666666664e-01f ) * t + 1.0f );
            c = d * ( ( ( ( ( -2.605e-07f * t + 2.47609e-05f ) * t - 1.3888397e-03f ) * t + 4.16666418e-02f ) * t - 4.999999963e-01f ) * t + 1.0f );
        }

        float math::tan( float a )
        {
            return tanf( a );
        }

        float math::tan16( float a )
        {
            float s;
            bool reciprocal;

            if ( ( a < 0.0f ) || ( a >= PI ) ) {
                a -= floorf( a / PI ) * PI;
            }
#if 1
            if ( a < HALF_PI ) {
                if ( a > ONEFOURTH_PI ) {
                    a = HALF_PI - a;
                    reciprocal = true;
                } else {
                    reciprocal = false;
                }
            } else {
                if ( a > HALF_PI + ONEFOURTH_PI ) {
                    a = a - PI;
                    reciprocal = false;
                } else {
                    a = HALF_PI - a;
                    reciprocal = true;
                }
            }
#else
            a = HALF_PI - a;
            if ( fabs( a ) >= ONEFOURTH_PI ) {
                a = ( ( a < 0.0f ) ? -HALF_PI : HALF_PI ) - a;
                reciprocal = false;
            } else {
                reciprocal = true;
            }
#endif
            s = a * a;
            s = a * ( ( ( ( ( ( 9.5168091e-03f * s + 2.900525e-03f ) * s + 2.45650893e-02f ) * s + 5.33740603e-02f ) * s + 1.333923995e-01f ) * s + 3.333314036e-01f ) * s + 1.0f );
            if ( reciprocal ) {
                return 1.0f / s;
            } else {
                return s;
            }
        }

        double math::tan64( float a )
        {
            return tan( a );
        }

        float math::asin( float a )
        {
            if ( a <= -1.0f ) {
                return -HALF_PI;
            }
            if ( a >= 1.0f ) {
                return HALF_PI;
            }
            return asinf( a );
        }

        float math::asin16( float a )
        {
            if ( FLOATSIGNBITSET( a ) ) {
                if ( a <= -1.0f ) {
                    return -HALF_PI;
                }
                a = fabs( a );
                return ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a ) - HALF_PI;
            } else {
                if ( a >= 1.0f ) {
                    return HALF_PI;
                }
                return HALF_PI - ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
            }
        }

        double math::asin64( float a )
        {
            if ( a <= -1.0f ) {
                return -HALF_PI;
            }
            if ( a >= 1.0f ) {
                return HALF_PI;
            }
            return asin( a );
        }

        float math::acos( float a )
        {
            if ( a <= -1.0f ) {
                return PI;
            }
            if ( a >= 1.0f ) {
                return 0.0f;
            }
            return acosf( a );
        }

        float math::acos16( float a )
        {
            if ( FLOATSIGNBITSET( a ) ) {
                if ( a <= -1.0f ) {
                    return PI;
                }
                a = fabs( a );
                return PI - ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
            } else {
                if ( a >= 1.0f ) {
                    return 0.0f;
                }
                return ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
            }
        }

        double math::acos64( float a )
        {
            if ( a <= -1.0f ) {
                return PI;
            }
            if ( a >= 1.0f ) {
                return 0.0f;
            }
            return acos( a );
        }

        float math::atan( float a )
        {
            return atanf( a );
        }

        float math::atan( float y, float x )
        {
            return atan2f( y, x );
        }
        float math::atan16( float a )
        {
            float s;

            if ( fabs( a ) > 1.0f ) {
                a = 1.0f / a;
                s = a * a;
                s = - ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
                    * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
                if ( FLOATSIGNBITSET( a ) ) {
                    return s - HALF_PI;
                } else {
                    return s + HALF_PI;
                }
            } else {
                s = a * a;
                return ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
                    * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
            }
        }

        float math::atan16( float y, float x )
        {
            float a, s;

            if ( fabs( y ) > fabs( x ) ) {
                a = x / y;
                s = a * a;
                s = - ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
                    * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
                if ( FLOATSIGNBITSET( a ) ) {
                    return s - HALF_PI;
                } else {
                    return s + HALF_PI;
                }
            } else {
                a = y / x;
                s = a * a;
                return ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
                    * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
            }
        }
        double math::atan64( float a )
        {
            return atan( a );
        }

        double math::atan64( float y, float x )
        {
            return atan2( y, x );
        }

        float math::pow( float x, float y )
        {
            return powf( x, y );
        }

        double math::pow64( float x, float y )
        {
            return powf( x, y );
        }

        float math::exp( float f )
        {
            return expf( f );
        }

        double math::exp64( float f )
        {
            return exp( f );
        }

        float math::log( float f )
        {
            return logf( f );
        }

        double math::log64( float f )
        {
            return log( f );
        }

        int math::bits_for_float( float f )
        {
            return static_cast< int >( fabs( static_cast< float >( ilog2( f ) ) ) ) + 1;
        }

        int math::bits_for_float( float max, float min )
        {
            int bitsForMax = static_cast< int >( fabs( static_cast< float >( ilog2( max ) ) ) ) + 1;
            int bitsForMin = static_cast< int >( fabs( static_cast< float >( ilog2( min ) ) ) ) + 1;

            return bitsForMax > bitsForMin ? bitsForMax : bitsForMin;
        }

        float math::fabs( float f )
        {
#if 1
            return fabsf( f );
#else
            int tmp = *reinterpret_cast<int *>( &f );
            tmp &= 0x7FFFFFFF;
            return *reinterpret_cast<float *>( &tmp );
#endif
        }

        float math::floor( float f )
        {
            return floorf( f );
        }

        float math::ceil( float f )
        {
            return ceilf( f );
        }

        float math::rint( float f )
        {
            return floorf( f + 0.5f );
        }

        float math::angle_normalize_360( float angle )
        {
            if ( ( angle >= 360.0f ) || ( angle < 0.0f ) ) {
                angle -= floor( angle / 360.0f ) * 360.0f;
            }
            return angle;
        }

        float math::fmod( float x, float y )
        {
            return fmodf( x, y );
        }

        float math::modf( float f, float& i )
        {
            return modff( f, &i );
        }

        float math::modf( float f )
        {
            float i;
            return modff( f, &i );
        }
    }   // ::math

}   // ::inf
